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IMPROVED AUTOMATIC COMPUTATION 
OF HESSIAN MATRIX SPECTRAL BOUNDS 

MORITZ SCHULZE DARUP AND MARTIN MONNIGMANN 


Abstract. This paper presents a fast and powerful method for the computation of eigenvalue 
bounds for Hessian matrices of nonlinear functions ip hi ^ —>• R on hyperrectangles 

B d lA. The method is based on a recently proposed procedure [9] for an efficient computation 
of spectral bounds using extended codelists. Both the approach from [9] and the one presented 
here substantially differ from established methods in that they do deliberately not use any interval 
matrices and thus result in a favorable numerical complexity of order 0{n) where N{ipi) denotes 

the number of operations needed to evaluate at a point in its domain. 

We improve the method presented in by exploiting sparsity, which naturally arises in the 
underlying codelists. The new method provides bounds that are as good as, or better than those 
from the most accurate existing method in about 82 % of the test cases. 

1. Introduction. We present important improvements for a recently proposed 
method (see 0) for the efficient calculation of spectral bounds for Hessian matrices 
on hyperrectangles. The improvements build on a systematic treatment of sparsity of 
the involved matrices, which will be shown to result in significantly tighter eigenvalue 
bounds. The problem can concisely be summarized as follows. Let : W C R" —>■ R. 
be a twice continuously differentiable function on an open set U C R" and let B = 
[xi,xi] X ■ • ■ X [Xj^,Xn] be a closed hyperrectangle in U. We seek bounds A, A C R 
such that the relations A < A < A hold for all eigenvalues A of all matrices H G 
{'S/^(p{x) I X G B}. More precisely, we solve the following problem: 

Find A, A G R such that A < min Amin(V^:/3(a:)) and max Ainax(V^(^(a:)) < A, (1.1) 

x^B x^B 

where XminiH) and Aniax(^f) denote the smallest and largest eigenvalue, respectively, 
of the symmetric matrix H G R"^". A bound A (resp. A) is called tight if there exists 
at least one x G B such that A = Ainin(V^vj(a;)) (resp. A = Xmax{X/'^<f{x))). Note that 
the problem statement (HI) does not necessarily imply that A and A are tight. 

Eigenvalue bounds A, A are used, for example, in numerical optimization methods 
to detect convexity, or to construct convex underestimators of nonconvex functions [U 
Ells]. If (11.111 yields A > 0 then (p is convex on the interior of the hyperrectangle B. 
While no conclusion on the convexity can be drawn if (HH) results in A < 0, the bound 
A can still be used to construct a convex underestimator for (p on B. Specifically, 

>fi{x) = p{x) - X A V , (Xi - Xi) {Xi - Xi) (1.2) 

is convex, coincides with p at the vertices of B, and bounds tp from below everywhere 
else in B. Since a large fraction of the total computation time is spent on the calcu¬ 
lation of convex underestimators in global optimization methods [1] , fast methods for 
solving (HI) are of interest. We briefly note that (ED must also be solved in certain 
problems in automatic control and systems theory. An illustrative example is given 

in 0- 

Existing approaches to solving ED proceed in two steps: First, a symmetric 
interval matrix (also called interval Hessian) that contains all Hessians X/'^(p{x) on B 
is calculated: 

Find , H = G R"^" such that < {y'^ip{x)).. < Hij 

for every f, j G {!,..., n} and every x G B. 
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This task can efficiently be carried out by combining interval arithmetics (lA, see 
m, for example) and automatic differentiation (AD, see [HUT], for example). In the 
second step, spectral bounds can be found by solving the following problem, which is 
similar to, but different from (ini): 

Find A, A S K such that A < min XminiH) and max XmaAH) < A, (1.4) 
~ ~ Hen Hen 

where V. = {H £ H = H^} is the set of all symmetric 

matrices that respect the bounds and H. Various approaches exist to solving (d 
(see, e.g., [2] O 0 [H])- However, since {V'^ip{x)\x £ B} £ V., problem (11.41) is 
conservative compared to the original problem m- In fact it is the very point of the 
method introduced in [9] and refined in the present paper to avoid computing interval 
matrices of the form d when solving d in order to avoid this conservatism. 

We briefly summarize the computational complexity of the existing methods. The 
computation of the matrices i/ and H in (11.31) requires 0{n^) N{tf) (resp. 0{n) N{(p)) 
operations if the forward (resp. backward) mode of AD is used, where N{(p) denotes 
the number of operations needed to evaluate at a point in its domain |4] . After the 
interval Hessian has been calculated, solving dd requires between O(n^) operations 
for the interval variant of Gershgorin’s circle criterion [HIS] and 0(2" n^) operations 
for Hertz and Rohn’s method pna. The latter method is an important benchmark 
in that it yields tight spectral bounds for the matrix set 'H. Albeit the conservatism 
in d, Hertz and Rohn’s method therefore provides the best possible option to 
solve d via d and d- 

The total numerical effort of any approach that uses d and d corresponds to 
the sum of the efforts for calculating , H and solving d- Thus, the numerical effort 
for the established methods varies between 0{n) N{(p) + 0{n^) (backward mode AD 
combined with Gershgorin’s circle criterion) and 0{n'^) N{(p) + 0(2" n^) operations 
(forward mode AD combined with Hertz and Rohn’s method). The major advantage 
of the direct method presented in |9] is its low computational complexity, which was 
shown to be of order 0(n) N{(p). 

It is the purpose of this paper to improve the method introduced in 0 such that 
sparsity can exploited to find tighter eigenvalue bounds. The improvements do not 
increase the numerical effort compared to the original method in |9]. In fact, sparsity 
needs to be investigated once during the automatic generation of the extended codelist. 
The computations required to evaluate the codelist to obtain eigenvalue bounds on 
a specific hyperrectangle B are no more expensive than those for the non-sparse case 
treated in [^. While the computational effort remains the same, the improved method 
results in significantly tighter eigenvalue bounds than the original procedure from [9] . 
To show this, we investigate 1522 examples (taken from the COCONUT collection [14], 
see m for details) and compare the eigenvalue bounds resulting from the improved 
procedure to those obtained with the original one [5] and to bounds obtained with the 
interval Hessian m and Gershgorin’s circle criterion and Hertz and Rohn’s method. 

We summarize the major aspects of the direct method for the computation of 
eigenvalue bounds from [9] in Sect. [3] Our main result, the exploitation of sparsity 
for the improvement of the eigenvalue bounds from [9], is stated in Sect. [4] We analyze 
1522 numerical examples from [TH US] in Sect. [^. Conclusions are given in Sect. 0 

^Results were obtained with Jcodegen, a code generator available from the authors on request, or 
to be used online on www.rus.rub.de/software/jcodegen. Jcodegen generates ANSI-C code for the al¬ 
gorithm described in Pron. ITTfl for a given function ip. In particular sparsity is treated automatically. 
The specific hyperrectangle B is passed to the resulting code as a runtime parameter. 
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2. Notation and Preliminaries. We frequently use index sets J C A/”, where 
M := and where := {i S N | m < i < n}. The complement of an index set 

J is defined as := M\J. The cardinality of an index set J is denoted by \J\. 

It is convenient to state eigenvalue bounds as intervals (e.g. A € [A, A]). Intervals 
[a, o] C M with a <a are further abbreviated by [o] := [o,a] whenever appropriate. 
Interval equality [a, a] = [6, b] , is understood to mean a = b and a = b. We frequently 
carry out calculations on intervals with standard interval arithmetics (lA) rules. The 
required rules are summarized in Lem. \TJ\ and Tab. 12.Il 

Lemma 2.1 (basic interval operations [IH])- Let [a] = [a, a] and [6] = [b,b] be 
intervals in R. Let a € [a,a], b G [b,b], and c S K 6e arbitrary real numbers. Then, 
the relations in the second column of Tab. \2.1\ hold under the additional restrictions 
stated in the last eolumn. 


Table 2.1 

Basic interval arithmetic. 


no. 

operation / bounds 

definition 

restriction 

1 

a + b 

G [a] + [b] 


= [a + b,a + b] 



2 

a b 

G [a] [b] 


= [mm{ab, ab,ab,a 

b}, niax{a b, a b, a b, a b}] 


3 

l/a 

G 1/H 


= [1/a, 1/a] 


0 i [a] 

4 

a"* 

£ [ar 


= [a"‘,a"*] 

if a > 0 or m odd 






= [a™-, a™-] 

if a < 0 and m even 






= ]0,max|a™',a’"|] 

otherwise 


5 

Va 

e yR 


= [y/a, Va] 


a > 0 

6 

exp a 

G exp ([a]) 


= [exp(a), exp(a)] 



7 

ln(a) 

G ln([a]) 


= [ln(a),ln(a)] 


a > 0 

8 

a + c 

G [o] + c 


= [a + c, a + c] 



9 

ca 

G c [a] 


= [ca, ca] 

if c > 0 






= [ca, ca] 

otherwise 



A lower case letter surrounded by brackets may refer to a real interval [x] = 
[x,x\ C M (as introduced above) or a hyperrectangle [x] = [x^,Ti] x • • • x [x„,Xn] C 
R" (with n > 2). In the latter case, the interval operations listed in Tab. 12.11 are 
understood to apply to every component. For a hyperrectangle [x] C R" and a 
nonempty index set J Q N with the m elements ji < • ■ • < jm, the term [xj] refers 
to the hyperrectangle [xj] = [x^J x • • • x [xj^]. 

It is furthermore convenient to use null matrices in R’”^’", which we denote by 
Om,r, when dealing with sparsity. For the special cases m = 0 or r = 0 we obtain an 
empty matrix. Formally, the empty square matrix 0o,o has no eigenvalues. It proves 
useful to assign the eigenvalue bounds [A] = [0, 0] to it. Finally, the Cartesian unit 
vector along the k-th direction is denoted by Cfc G R". 


3. Direct computation of eigenvalue bounds for Hessian matrices on 
hyperrectangles. We summarize the method introduced in for the direct solution 
of (II.1|) as needed in the present paper. We assume the function (p can be evaluated 
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at an arbitrary point x £lAhy carrying out a finite sequence of operations of the form 


2/1 = 

Xl 


2/n = 

Xn 


Vn+l = 

yn+2 = 

^n+i{yi, ■ 
^n+2{yi, ■ 

, (3.1) 

■-,2/n,yn+l) 

Vn+t — 


■ ■ ; Vni ^n+l: ■ - ■ i ^n+t—l) 

p = 

Vn+t 


where each k = 1,... 

, t, represents one of the elementary operations listed 


in the first column of Tab. o We treat the same operations as in [5] for ease of 
comparison. Note that additional unary operations can be added according to the 
rules given in [5]. We refer to (13.ip as the codelist of the function ip. 

The codelist (ED can be used to evaluate the function value >p{x) at a specific 
point X in its domain. Using automatic differentiation (AD) [11] the codelist (ED 
can be extended in such a way that the gradient V(p{x) or the Hessian V'^(p{x) at the 
point X are calculated. Moreover, using AD and interval arithmetic (lA), (13.11) can 
be modified such that interval extensions, interval gradients or interval Hessians of p 
on hyperrectangles B <Z lA are computed. In fact, extended codelists are commonly 
used to solve problem ED as part of the established procedures for the computation 
of eigenvalue bounds (see, e.g., |2]). In contrast, the method introduced in only 
requires the interval gradient, but not the interval Hessian. Essentially, the codelist 
is extended by arithmetic operations that compute the eigenvalue bounds for the 
Hessian of the intermediate function in every codelist line. Formally, this leads to the 
extended codelist which we introduce in the following theorem. 

Theorem 3.1 (algorithm for direct eigenvalue bound computation jH] Prop. 
4.2]). Assume tp is twice continuously differentiable on U and can be written as a 
codelist ED- Let B = [xi] X • • • X [xn] dU be arbitrary. Then, for all x G B, we 
have p{x) G [p], Vp{x) G [Vt/s], and [Amin(V2(p(a;)), Amax(V2(/5(a;))] C [A,^], where [p], 
[Vt/?], and [Ay,] are calculated by the following algorithm. 

1. Fork = l,...,n, set [yk] = [xi,,Xk], [Vyk] = [efc,efc], and [Afe] = [0,0]. 

2. For k = n + 1,... ,n + t, calculate [yk], [Vj/fe] and [Afe] according to the third, 
fourth, and fifth column of Tab. \3.1\ 

3. Set [p] = [vn+t], [V(^] = [Vyn+t], and [Ay,] = [A„+t]. 

In Tab. 13.11 we use the interval operators [As([a])] and [At([a], [6])], which are 
defined according to 


[A.([a])] = 
and [At([a], [5])] = 


[a]^ if m = 1, 

[0 1 Sil 1 niax{a^, af }] otherwise, 

2 [a] [6] if m = 1, 

[-/3, 13] + { 3 L^,a^] [bi,bi] otherwise 


(3.2) 

(3.3) 


for hyperrectangles [a], [b] C M™, where /3 = (X]™ i maxja^, af}){jy^i xaa.x{lg , b^}). 

We refer to (SJ Lems. 2.2 and 2.3] for details on [As([a])] and [At([a], [6])]. 

4. Improved computation of eigenvalue bounds using sparsity. If spar¬ 
sity is exploited, tighter eigenvalue bounds can be obtained than those that result from 
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Table 3.1 

Rules for the calculation ofyi^, [yk], [Vj/fe] and [Aj.] in the k-th line of the codelist 113.Ill . [Vj/fe] 
refers to the interval gradient of line k with respect to x. The interval operators [As([a])] and 
[At([a], [6])] are defined in 113.2|| and II3.3II . 


op 

Vk 

[yk] 

[Nyk] 

[Afe] 

var 

Xk 

[Xk] 

\^k 1 

[0,0] 

add 

Vi + Vj 

[yi] + [yj] 

Nvi] + [V 2 /,] 

[Ai] + [Xj] 

mul 

yiVj 

[yi] [yj] 

[2/4][V2/i] + [2/i][V2/j] 

[yAM + [yi][>^j] + [At([V 2 /i], [Nyj])] 

powNat 

yT 

[yr 

m [yir-^ [Vyi] 

m[yir~'%{m-l)[As{[Vyi])] + [yi][Xi]) 

oneOver 

i^/yi 

i/[j/i] 

-[yk]^ [Nyi] 

h/kf (2 [yk] [As([V 2 /i])] - [Ai]) 

sqrt 

Vyi 

IVW]] 

1/(2 [yk]) [Nyi] 

1/(2 [2/fe])([A.] + l/(-2[2/,])[A,([V2/i])]) 

exp 

exp( 2 /i) 

[exp([ 2 /i])] 

[yk] [Nyi] 

[yk] ([AsQVi/i])] + [Ai]) 

In 

ln( 2 /i) 

[ln([ 2 /i])] 

l/[j/i] [V2/i] 

1/bi] ([I'i] - 1/bi] [As([V2/i])]) 

addC 

yi + c 

[yi]+[c,c] 

[V2/»] 

[Ai] 

mulByC 

cyi 

c[yi] 

c[V2/i] 

c[Ai] 


the method summarized in Section [3] This is evident from the following motivating 
example. 

Example 1 (method from [5] applied to (fi^) = xf + x^)- Consider the function 
—>■ R with (p{x) = xf + x^- Theorem \3.1\ results in the following extended 
codelist. Note that the expressions for yp listed in (14.11) do not result from Thm. \3.1\. 
but are only given for illustration of the codelist (EHI) of ip. 


k 

yk 

[yk] 

[Nyk] 

[Afe] 

1 

XI 

bi] 

[ei,ei] 

[ 0 , 0 ] 

2 

X2 

[X2\ 

[62,62] 

[ 0 , 0 ] 

3 

yl 

[yi? 

2 [2/1] [V2/1] 

2 ([As([V 2 /i])] + [2/1] [Ai]) 

4 

yl 

[y2? 

2 [2/2] [V2/2] 

2 ([As([V 2 / 2 ])] + [2/2] [A2]) 

5 

2/3 + 2/4 

[2/3] + [2/4] 

[Vi/s] + [V2/4] 

[As] + [A4] 


= 2/5 

[ t ] = [2/5] 

[Vvp] = [V2/5] 

[Aip] = [As] 


Evaluating the extended codelist (14.11) for the hyperrectangle B = [0,1] x [0,1] by 
computing [yp], and [Afc] and storing the results line by line yields 


H = [0,1], [Vt/i] = ([l,l],[0,0]r, [Ai] = [0,0], 

[7/2] = [0,1], [Vy2] = ([0,0],[l,l])^, [A2] = [0,0], 

[2/3] = [0,1], [Vj/ 3 ] = ([0,2l,[0,0ir, [A3] = [0,2], 

[2/4] = [ 0 , 1 ], [V2/4] = ([ 0 , 0 ], [ 0 , 2 ]f, [A4] = [ 0 , 2 ], 

[2/5] = [0,2], [V2/5] = ([0,2],[0,2]r, [As] = [0,4], 


(4.2) 


where [As([Vj/i])] = [0,1] and [As([V 2 / 2 ])] = [0,1] according to Eq. ([3)3]). Thus, we 
obtain the eigenvalue bounds [A,^] = [A 5 ] = [0,4] for 'V‘^{p{x)) on B. Now, consider 
the functions 5 , —)• R with g{x) = x\ and h{x) = a;|. From 


V^g{x) = ^ ^ H ^ and V^h{x) = ^ ^ , 


we infer that both V^g{x) and V^h{x) have the eigenvalues 0 and 2 for every x G B. 
Hence, the eigenvalue bounds [Ag] = [A 3 ] = [0,2] and [Ati] = [A 4 ] = [0,2] that result 
in line 3 and f of extended codelist (ED are tight. The eigenvalue bounds [A 5 ] = 
[A 3 ] + [A 4 ] = [0,4] that result in the subsequent line are conservative, however. In fact. 
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the Hessian of ip reads 


2 0 
0 2 




V^g{x) + V^h{x) = 


for all X £ B and the tight eigenvalue bounds obviously read [A*] = [2,2]. 

The Hessian matrices V^g(x) and \/^h{x) in Exmp.[T]have zero eigenvalues which 
disappear when adding the two functions to f(x) = g{x) + h(x). The situation 
illustrated in Example [T] arises naturally in the codelists introduced in Section [31 
because codelists build up functions of many variables from functions of very few of 
these variables. In order to mitigate eigenvalue bound overestimation in these cases, 
we need to consider functions like g{x) and h{x) in Exmp. [T] as functions of only 
those variables that they actually depend on nonlinearly. To this end, some simple 
terminology and intermediate results are introduced in Sect. 14.ll Subsequently, sparse 
sums, products, and compositions are treated in Sects. 13131 14.31 and 14.41 respectively. 
Section 14.51 summarizes how to compute the improved eigenvalue bounds based on 
the rules introduced in Sects. I4.2H4.41 

4.1. Sparsity handling using reduced Hessians and reduced gradients. 

As pointed out in Exmp. Ill sparsity occurs if functions depend at most linearly on 
some variables Xi, where at most linear dependence is defined as follows. 

Definition 4.1 (at most linear dependence). Let f : U ^ M. be a continuously 
differentiable function on an open set li C R". Let i £ J\f. The function f is said to 
depend at most linearly on Xi if there exists a c S R such that 

df 

— — (x) = c for all X £14. (4-3) 

UXi 

The function f obviously is independent of Xi if (lOI) holds with c = 0. 

Assume a function / is known to depend at most linearly on Xj for all j £ Cf, 
where Cf C A/” is a given index set. Then, only the eigenvalues of the reduced Hessian 
(see Def. 14.21) associated with the index set J = C^j are nontrivial, i.e., not necessarily 
equal to zero. 

Definition 4.2 (reduced Hessian V^/(a;)). Let f : 14 ^ M. be a twice con¬ 
tinuously differentiable function on an open set 14 C R". Let J £f M he an index 
set and let m = \J\. If m = 0 set V^/(a:) = Oq.o, otherwise denote the m ele¬ 
ments of If by ji < • ■ ■ < jn-m in ascending order and define the reduced Hessian 
V^jf{x) £ })y iig elements 




d'^fjx) 
dxj.dxj^ ’ 


where i,k £ Ni^m- 

We also need to consider reduced gradient vectors. 

Definition 4.3 (reduced gradient Vy/(x)). Let f : U ^ R be a continuously 
differentiable function on an open set 14 C R". Let fl <£ M be a nonempty index set 
and let m = \J\. Denote the m elements of J by ji <■■■ < jn-m in ascending order 
and define the reduced gradient V jf{x) £ R"“”^ by its elements 


where i £ 


(Vy/(x)), 


dfjx) 

dxj, ’ 
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Note that J may be empty in Def. 14.21 while there must exist at least one element 
in J in Def. 14.31 This difference arises since codelist lines may depend at most linearly 
on all variables Xi but they are never independent of all Xi. 

We can easily evaluate eigenvalue bounds for the Hessian of a function from 
eigenvalue bounds for its reduced Hessian. This is stated precisely in Lem. 14.41 

Lemma 4.4 (spectral bounds for Hessian from reduced Hessian). Let f denote a 
twice continuously differentiable function f :U on an open set U. Let the index 
set Cf C J\f be such that f depends at most linearly on Xi for all i € Cj. Let B CU 
and let the interval [Ay] C M. be such that 


AI < minAi„in(V£c/(x)) and maxAi„ax(V£c/(x)) < x\. 
•’ x^B J x£B f 

Then, the eigenvalues of the Hessian V^/(x) on B lie in the interval 


(4.4) 


[A^ 


/J 


[A/]= [0,0] 


*/£y = 0 , 

ifCf=M, 


(4.5) 


[min{Ay, 0}, max{Ay, 0}] otherwise. 


Proof We consider the cases in (14.51) separately. Cf = % implies C‘j = Jf and 
consequently V^/(x) = V‘j^cf{x), which proves the first case. In the second case, i.e., 
Cf = Af, we have V^/(x) = S7j^^f{x). Since / depends at most linearly on Xi for 
all i G Cf, we find V£^/(x) = 0n,n- Thus, the eigenvalue bounds [A/] = [0,0] hold. 
Regarding the third case, we note that % C Cf C N implies 0 C £y C A/”. Thus, 
m = \Cf \ satisfies 0 < m < n. Without loss of generality we assume Cf = Ni,m- 
Then 


Now consider an arbitrary but fixed x G B. We obtain 

Amin(V^/(x)) = min{Ainin(V£c/(x)),0} and (4.6) 

Amax(VV(a^)) = max{An,ax(V|c/(x)),0} (4.7) 

based on the block-diagonal structure of V^c/(x). Bounding (|4.6p below and bound¬ 
ing (14.71) above for all x G H yields 

[minAmin(V^/(x)),maxAi„ax(V^/(x))] C [min{At, 0}, max{At, 0}] 
xSB xeB ■> ■’ 

according to Eqs. (ITBI) and (14.7p and condition (14.4L □ 

4.2. Improved eigenvalue bounds for the sum of two functions. We col¬ 
lect some recurring conditions first. 

Conditions 4.5. Let g and h denote twice continuously differentiable functions 
g : U ^ M. and h '.lA ^ on an open set lA C K”. Let the index sets Cg C J\f and 
Ch C J\f be such that g (resp. h) depends at most linearly on Xi for all i G Cg (resp. all 
i G Ch). Moreover, let the index sets Ig C Cg and Ih Q Ch with Ig C Af and Ih G N 


V^fix) = 


S 


0 „- 


n—m.m 


vL/(x) 


0m,n—71 


On- 


n—m,m 


On—m,n—771 
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be such that g (resp. h) is independent of Xi for all i € Ig (resp. all i G I}^). Let 
B <ZU and assume there exist intervals [Aj] C K and [Aj^] C K. such that 


AJ < minAmin(V|cg(a;)) 

A|i a niin Amin(^£c/?-( x)) 
x^B ^ 


and maxAinax(V|cg(x)) < A^, 

xGB b b 

and max Aniax(V£c/i(a;)) < aJ^. 


(4.8) 

(4.9) 


Now assume Conds. 14.51 hold and we intend to calculate eigenvalue bounds for 
V^/(x) on a hyperrectangle B for f{x) = g{x)-\-h{x). We could determine eigenvalue 
bounds for the full Hessians \/^g{x) and V‘^h{x) with Lem. 14.41 and apply the rule for 
the eigenvalue bounds of the sum of full Hessians (line add in Tab. 13.11 reproduced 
from i)- However, we show in Lem. 14.81 below that it is advantageous to, roughly 
speaking, carry out calculations with the sparse Hessians as long as possible and to 
apply Lem. 14.41 as late as possible. We first state the rules for determining jC,f,If and 
the eigenvalues of the reduced Hessian of / in Lems. 14.61 and 14.71 respectively. The 
trivial proof of Lemma 14.61 is omitted for brevity. 

Lemma 4.6 (index sets for sums). Assume Conds. \4.5\ hold and consider the 
function f : U ^ M. with f{x) = g{x) + h{x). Let Cf = Cg A Ch and Xf = Xg nXh- 
Then, f depends at most linearly on Xi for all i G Cf and f is independent of Xi for 
all i G Xf. 

Lemma 4.7 (spectral bounds for reduced Hessian of sums). Assume Conds. \4-5\ 
hold and consider the function / : —>■ K, f{x) = g{x) + h{x). Let Cf — CgG\ Ch. 
Then, 

At < minAi„in(V£c /(a;)) and max Ai„ax(V£c/(a;)) < aI, (4.10) 

x^B f x^B f 

where [Aj] is computed according to the rules listed in Tab. \4-l\ 


Rules for the computation of eigenvalue bounds [aJ.] for the reduced Hessian V^c /(jr) of a sum 

f{x) = g(x) + h(x). Let Cu be short for Co ■= Cg U C},. See the end of Sect. f ..2 for a discussion 
of the eight cases. 


case 

[At] 

condition 

1 

[0,0] 

Cg=N ACh=N 

2 

[At] 

Cg G N A Ch = M 

3 

[At] 

Cg=N AChGN 

4 

[miniAt, A),}, max{ At, aI}] 

Cg G N A Ch G N A Cu = N 

5 

[At] + [At] 

Cu GM ACg = Ch 

6 

[At] + [mm{At, 0 }, max{Al, 0 }] 

Cu G Af A Cg G Ch 

7 

[min{At, 0 }, max{At, 0 }] + [A],] 

Cu G N A Ch G Cg 

8 

[minlAt, 0 }, max{At, 0 }] + [min{At, 0 }, max{Al, 0 }] 

CoGAT A Cg<^Ch A Ch^Cg 


Proof. We prove the fourth case in Tab. 14.11 since it will be instrumental for 
Exmp.[^ All other cases in Tab. I4.1| can be proven analogously. The reduced Hessian 
of / reads Vj.cf{x) = yj^ag{x) + Vj.ch{x). From Cg G N, Ch G N, and Cyj = 

CgG Ch = M, we infer tf) g Cg G N and 0 C C Af and consequently 0 C C A/” 
and % G Cf G N. Thus, the cardinalities r = \Cg\ and s = \Cf\ satisfy 0 < r < n 
and 0 < s < n. Moreover, Cg U Ch = Af implies CfCi Cf = 0. Hence, there does not 






exist any index i G J\f such that both * S and i G We assume = Ni_r and 
= Nr+i,r-+s without loss of generality. Note that implies = Ni_f.+s 

and m = |£'i| = r + s under this assumption. Thus, V£o/(a;) equals 

J 


( V|eff(a;) 0r,s \ 

\ 0s,r. 0s,s / V / y 0s,T- '^C'^h{x) j 

The block-diagonal structure implies 

'^min(^ £■= /(^)) ~ ^ill{'^min(^'^min /l(x))} and 

f g h 

'^max (^f ix')') = max! Aniax(^ C^9^x') \ Amax(^^(^))} 

/ 9 h 


(4.11) 


(4.12) 

(4.13) 


for an arbitrary but fixed x G B. Bounding (14.121) below and bounding (14.131) above 
for all X S S yields 


[min Amin( f{x)), max Amax 

x^B x€B 


(VV(a;))] C [min{A^,At},max{Ag,Al}] 


where we used Eqs. (14.121) and (14.131) and Conds. 14.51 Thus, the eigenvalues of 
V£<=/(x) on B lie in the interval [Aj] = [min{Ag, Aj^}, max{Ag, Ajj}] as claimed in 
Tab. HU □ 

We anticipated the bounds from Lem. 14.71 can be shown to be as tight as or 
tighter than those from the original method proposed in [5] that does not account for 
sparsity. This can now be shown in Lemma H751 below. Recall the bounds in [S] result 
in [A/] = [Ag] + [A/i] for f{x) = g{x) + h{x) according to O Prop 3.2.(iii)]. 

Lemma 4.8 (improved bounds for sums). Assume Conds. \4-5\ hold and let f, 
[Aj], and Cj be as in Lem.\4^ Let [A/], [Ag], and [A^,] be the eigenvalue bounds for 
the Hessians V^/(x), 'V^g{x), and 'V^h{x) on B, calculated according to Eq. (|4.5I) . 
Then, 


[A/] C [Ag] + [A?i]. 


(4.14) 


Proof. We prove the relation for the fourth case in Tab. 14.11 The remaining 
cases can be proven analogously. As pointed out in the proof of Lem. 14.71 we have 
0 C £g C A/" and 0 C £/» C A/". Thus, the r.h.s. in (14.141) yields 

[Ag] + [Xh] = [min{Aj,0},max{Ag,0}] + [min{Aj^, 0}, max{Al, 0}] 

= [min{Ag, 0} + min{A|j, 0}, max{Ag, 0} + maxjAjj, 0}], 

= [min{Aj + A|^,A],,A|^,0},max{Ag + aI, Ag, Ai,0}], (4.15) 

where the equations hold according to the third case in (14.51) . by definition of the sum 
of two intervals (see Tab. \2A\i . and by definition of min{-} and max{-}, respectively. 
To evaluate the l.h.s. in (I4.14L we have to analyze the index set £/. We obviously 
have £/ = CgfiCh C Af. Thus, the second case in Eq. (14.51) does not apply. However, 
from the conditions characterizing the fourth case in Tab. 14.11 it is not clear whether 
£/ = 0 or D 0. Thus, according to (14.51) . the l.h.s. in (14.14p results in 


[A/] = 


[A}] _ if £/ = 0, 

[min{Aj, 0}, max{Aj, 0}] if0Cl3/CA/’. 
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However, since [a] C [min{a, 0}, max{a, 0}], the relation [A/] C [min{A^, 0}, max{Aj, 0}] 

holds in both cases. Since [Ajp] = [niin{Ag, Ajj},max{Ag, aJ^}] according to Lem. H771 
(resp. Tab. 14.11) . we obtain 

[A/] C [min{min{A^,Aj;},0},max{max{Ag, a1},0}], 

= [min{A^,Aj;,0},max{Ag, A^O}]. (4.16) 

Comparing Eqs. (14.151) and (14.161) yields 

[min{Aj,Aj),0},max{Ag, a 1,0}] C [min{Aj + Aj^, A^, A],, 0}, max{Ag + aI, a],, a],, 0}], 
which proves (14.141) . □ 


Table 4.2 

Fundamental cases underlying Tab. \4.1\ 


case 

reduced Hessian 

contribution 

condition 

(i) 

o 

o 

o 

none 

£= = 0 A £); = 0 

(ii) 


first Hessian 

£^ D 0 A £^ = 0 

(iii) 

Vlgh{x) 

second Hessian 

£g = 0 A £^ D 0 

(iv) 

Vlcg{x) + Vlghix) 

both Hessians 

£^ D 0 A £^ D 0 


Lemmas 14.71 and 14.81 are based on the eight cases listed in Tab. 14.11 Since they 
are not obvious at first sight, it is instructive to see how these eight cases arise from 
the four simpler ones listed in Tab. 14.21 The first case in Tab. 14.21 applies if both 
V‘j^cg{x) and \/j^ch{x) vanish because of Cg = C'^ = 0. Since these two conditions, 
i.e. = 0 A £)( = 0, are equivalent to £g — Af A Ch = Af, case (i) in Tab. 14.21 is 
equivalent to case 1 in Tab. 14.11 Analogously, cases (ii) and (iii) in Tab. 14.21 where 
either V'jrcg{x) or \7‘j-ch{x) contribute to V'jrcf{x), are equivalent to cases 2 and 3 in 
Tab. 14.11 respectively. It remains to relate case (iv) in Tab. l4.21 to cases 4-8 in Tab. 14.11 
In fact, the conditions of the cases 4-8 in Tab. 14.II all imply D 0 and £((3 0, which 
are the defining conditions for case (iv) in Tab. 14.21 Figure HIT] illustrates that every 
instance of case (iv) from Tab. 14.21 actually uniquely belongs to one of the cases 4-8 
from Tab. 14.11 


J 


□ 




case 4 


case 5 



□ 


case 6 


J 



case 7 



case 8 


Fig. 4.1. Illustration of the sparsity patterns of the reduced Hessians and ^'j^ch{x) 

for the cases f-8 in Tab. \4. and Lem. [13 White areas correspond to zero blocks in the Hessians, 
black areas to non-trivial blocks. The black areas in case 4 correspond to the non-trivial submatrices 
and S7\ah(x) on the l.h.s. of Eq. 114.till . 


4.3. Improved eigenvalue bounds for the composition of two functions. 

We collect some recurring conditions again first. 

Conditions 4.9. Assume Conds.^^hold. Letr : V —>■ M 6 e a twice differentiable 
function on an open set V C M and assume g{x) G V for every x & Li. Moreover, 
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assume there exist intervals [r'\ C M and [r"] C M and a hyperrectangle [Vg] C M" 
such that 


r' < r'{g{x)) < r', r" < r"{g{x)) < r", and V 5 . < (Vg(x))i < Vg, (4.17) 

for every x € B and every i € M, where r'{z) and r"{z) refer to the first and the 
second derivative ofr{z), respectively. 

The following lemma, which we state without proof, provides rules for the iden¬ 
tification of at most linear dependencies and independencies of compositions. 

Lemma 4.10 (index sets for compositions). Assume Conds. [4-9\ hold and con¬ 
sider the function f : U ^ M., f{x) = r{g{x)). Let 




and 


j Cg if r is an affine function, 
\ Xg otherwise. 


(4.18) 


Then, f depends at most linearly on Xi for all i G Cf and f is independent of Xi for 
all i G If. 

Bounds for compositions can now be calculated as follows. 

Lemma 4.11 (spectral bounds for reduced Hessian of compositions). Assume 
Conds. \4.9\ hold and consider the function / : W —>■ K, f{x) = r{g{x)). Let Cf be 
defined as in Lem. \4.10\ Then the bounds (I4.10p hold for [Aj] computed according to 
the rules listed in Tab. \4.3[ 


Table 4.3 

Rules for the computation of eigenvalue bounds [aJ.] for the reduced Hessian V^c fO of com¬ 
positions f{x) = r{g{x)). 


case 

[x\] 

condition 

1 

n [a.([V£=5 ])] 

Cg=N 

2 

n [W([V£.5])] + [r'] [At] 

Cg CAf ACf =Cg 

3 

[r"] [As([V£=g])] -b [r'] [min{At, 0}, max{Ag, 0}] 

Cg GN /\Cf C Cg 


Proof. We prove the last case in Tab. 14.31 The remaining cases can be proven 
analogously. The reduced Hessian of / reads V^c/(a:) = r"{g{x)) V£= g(a;) V'^cg(x) -b 

r'{g{x)) Vj-c g{x). Combining the two conditions of case 3 in Tab. l4.3l vields Cr CCgC 
M, which implies C'f D £g D 0. Thus, m = |£j| and r = |£g| satisfy m > r > 0. We 
assume and Cg = without loss of generality, tinder these assumptions, 

we obtain 


= r"{g{x)) Vcp{x) Vj. c/(x) + r'{g{x)) 


^ Cgdi^') \ 

ds,?" ds,s / 


(4.19) 


where s = m — r > 0. Since 'Vg{x) G [Vg] for every x G B, we find 

As([V£=5 ]) < C}g{x) Vjcg(x)) and Xmux{X7cpix) vjcg(x)) < As([V£cg]), 

for every x G B according to [SI Lem. 2.2]. Combining this intermediate result with 
the bounds on r'{g{x)) and r"{g{x)) from Conds. [DU yields the eigenvalue bounds 

[X\] = [r"] [A«([V£=g])] + [r'j [min{A^, 0}, max{Ag, 0}] 
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on B. □ 

Lemma r4.12l below shows that the bounds from Lem. 14 .11] are as tight as or tighter 
than those from the original method proposed in [^. Recall the bounds in result 
in [A/] = [r"] [As([Vg])] + [r'] [Ag] for f{x) = r{g{x)) according to O Prop 3.4]. 

Lemma 4.12 (improved bounds for compositions). Assume Conds. hold and 
let f, [Aj], and Cf be as in Lem. Let [A/] and [Ag] he the eigenvalue bounds for 

the Hessians V^/(a;) and V‘^g{x) on B, calculated according to Eg. (14.51) . Then, 

[A/] C [r"] [A,([Vg])] + [r'j [Agj. 


Since the proof is very similar to the proof of Lem. 14.81 we omit it. 

4.4. Improved eigenvalue bounds for the product of two functions. We 

begin by collecting recurring conditions again. 

Conditions 4.13. Assume Conds. |.^.5| hold and assume there exist intervals 
[g] C M and [h] C M and hyperrectangles [Vg] C M" and [V/i] C M" such that 

9 < 9{x) <g, h< h{x) < h, Vg^ < {Vg{x))i < and Vfe .- < (V/i(a:))i < V/ii 

for every x G B and every i G Af. 

The following lemma provides rules for the identification of at most linear depen¬ 
dencies and independences of products. 

Lemma 4.14 (index sets for products). Assume Conds. \4-18\ hold and consider 
the function f : U ^ M., f{x) = g{x) h{x). Let If = Ig Dlh and Cf = Ig nlh- Then, 
f depends at most linearly on Xi for all i G Cf and f is independent of Xi for all 

*el/. 

Based on Conds. 14.131 and Lem. 14.141 we are able to compute bounds on the 
spectrum of f{x) according to the rules summarized in Lem. 14.151 and Tab. 14.41 
As a preparation, we introduce the interval operators 


[Ar([a], [ 6 ])] := [min{a, 6 }, max{a, &}] 
1 


and 


[A*([a],[ 6 ],[c])] := 


-I- 6 — \/{a — b)^ + d, a + b+ {a — b)"^ + d 


(4.20) 

(4.21) 


for real intervals [a], [6], [c] C M, where d = 4 max{c^, c^}. Definition (14.201) is only 
introduced for the sake of a compact notation. Whenever it is more instructive, we 
use the notation on the r.h.s. of (I4.20|) . 

Lemma 4.15 (spectral bounds for reduced Hessian of products). Assume Conds. \ 4 - 13 \ 
hold and consider the function / : W —>■ M with f{x) = g(x) h(x). Let Cf be defined as 
in Lem. \4-14\ Then the bounds (14.101) hold for [A^j computed according to the rules 
listed in Tab. \4-4\ 

Proof. We prove case 10 from Tab. 14.41 Cases 4 and 7 can be shown analogously. 

The remaining cases can be proven in the same fashion as those treated in the proofs 
of Lems. l4Tl and l4.lT] The reduced Hessian of /, which reads 

V£=/(a;) = V£jg(a;)Vj=/i(a;) -b Vc=/i(a;)Vjcg(a;) -b h{x) V^cg^x) (4.22) 

in all cases, is a two-by-two matrix with a particularly simple block structure in 
case 10. To see this, first note that g and h are independent of all but one variable 
each (the conditions \Ig\ = n — 1 and \Ih\ = n — 1 imply IX^j = 1 and \If \ = 1 ). 
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Table 4.4 

Rules for the computation of eigenvalue bounds [aI] for the reduced Hessian 'V^cf{x) of a 




product f{x) = g{x)h{x). The expressions [At], [A^^oj? [■^h.o]? o,re short for [At] = 

[At([V£=g], [V£=h])], [Ag,o] = [min{At,0},max{Ag,0}], [A^^q] = [min{At, 0}, max{Aj,, 0}], Cu = 
Cg U Ch and C<-\ = Cg Pi Ch- Condition C* reads {Xg LlXf^ = Af) A {\Xg\ = n — 1) A = n — 1). 


1 [At] 

2 [At] + [/i] [A^] 

3 [At] + [/i] [A9 ,o] 

4 [A.([/t][At],[0,0],[Vi=5][Vi.fe])] 

5 [At] + [<?] [Al] 

6 [At] + [<;] [A/i,o] 

7 [A*([0,0],M[At],[Vxcg][Vicfe])] 

8 [At] + [A,([A.([fe][At],[5][Al])],[0,0])] 

9 [At] + [A.([/t][At],[s][At])] 

10 [A*([/t][At],[5][Al],[Vx=<7][Vxj^])] 

11 [At] + [h] [At] + [ff] [At] 

12 [At] + [/i] [A^] + [g] [Ah.o] 

13 [At] + [/i] [A 9 ,o] + [g] [Aj,] 

14 [At] + [A.([/t][At] + [5][At].[0,0])] 

15 [At] + [Ar([/i] [A^j + [g] [Ah.o], [0, 0])] 

16 [At] + [Ar([/i] [Ag,o] + [g] [A]j], [0,0])] 

17 [At] + [/i] [Ag.o] + [g] [Ati,o] 


condition 


Cg=N /\Ch=M 

Cg CAT ACh ACf ^ Cg 

Cg cAf A Ch =Af A Cf cCgA-iCi, 

Cg CAT A Ch= AT A C^, 

Cg=Ar A ChCAf A Cf=Ch 
Cg = Af A Ch C Af A Cf C Ch A -iCic 
Cg=Af A ChCAf A Ci, 

Cg C Af A Ch C Af A Cu = Af A Cf C Cn 
Cg <Z Af A Ch <Z Af A Cu = Af A Cf = £n A-.(7* 
CgCN AChCN AC^, 

Cu G Af A Cf = Cg = Ch 
Cu cAf ACf = Cg C Ch 
CuGAf ACf = ChG Cg 
CuGAf ACf GCg^ Ch 
CuGAf ACf G CgG Ch 


CuGN ACg(^ChACh<^Cg 


Moreover, IgUlh = Af implies Ig = 0, which implies g and h depend on two 
different variables. Without loss of generality we assume g depends on xi, and h 
depends on X 2 , i.e., 1° = {1} and = {2}. As a further preparation note that 
£/ = Ig n Ih, which holds according to Lem. 14.141 implies Cj = Ig L) which 
evaluates to = {1,2}. Since = {1,2} and g only depends on xi (resp. h only 
depends on X 2 ), we have 


'^Cpix) = 

respectively 

VctKx) = 


£:9{x)\ 

£^9ix)J 






^%9{x) = 




0 


Substituting (14.231) and (14.241) into (14.221) yields 


^%f{x) = 





c-fHx) = 


hMx) 


dx-i 


9{x)^h{x 




9{x)^M^) 


(4.23) 


(4.24) 


where all entries are scalars. Now, consider the matrix set 

n = {HG k2x2 I^ ^ 


and observe {V£c/(x) G G B} G %. To see this, note that {V|cg(a;) G 

M[x G B] C [AJ] and {V£c/i(a;) G M | x G S} C [Aj^], since the eigenvalue of a 
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matrix M € is A = Mi i. According to Lem. lA.ll stated in the appendix, 

eigenvalue bounds for the matrix set T-L and consequently for V£c/(a;) on B read 

[A*([/i] [Aj], [g] [Ajj], [Vic 5 ][Vi=/i])] as claimed in Tab. S^l 0 

Lemma [4. 161 shows that the bounds from Lem. 14.131 are as tight as or tighter than 
those from the original method proposed in j9]. Recall the bounds in result in 
[A/] = [At([V 5 ], [V/i])] + [h] [Ag] + [g] [A/^] for f{x) = g{x) h{x) according to (9] Prop 
3.2.(iv)]. We omit the proof of Lem. 14.161 since it is similar to its counterparts in 
Sect. 14.21 

Lemma 4.16 (improved bounds for products). Assume Conds. \4-9\ hold and let 
f, [A|], and Cj be as in Lem. \4-lT\ Let [A/], [Ag], and [A?i] be the eigenvalue bounds 
for the Hessians f{x), 'S/^g{x), and V‘^h{x) on B, calculated according to Eg. (|4.5I) . 
Then, 


[A/] C[A*([Vg],[Vh])] + [h] [Ag] + [ 5 ] [A^]. 


4.5. Numerical computation of improved eigenvalue bounds. In this sec¬ 
tion, we combine the results from Sects. HTT] through l4.4l in order to compute improved 
eigenvalue bounds using a codelist. Formally, this leads to the extended codelist in 
Prop. 14.171 

Proposition 4.17 (algorithm for the computation of eigenvalue bounds using 
sparsity). Assume ip is twice continuously differentiable on U and can be written as 
a codelist iiD with t €N operations. Let B = [xi] x • • • x [xn] dU be arbitrary. Then, 
for all X d B, we have ip{x) G [p], 'Vp{x) G [Vt/s], and [A„iin (V 2 <^(a;)),A„,ax(VV(a:))] C 
[Ag,], where [t/s], [Vt/?], and [A,^] are calculated by the following algorithm. 

1. Fork = l,...,n, setik =M\{k}, Ck = Af, [yk] = [Xk,Xk], ['^yk] = [ek,ek], 
and [A^] = [ 0 , 0 ]. 

2. For k = n + 1,... ,n + t, evaluate Ik and Ck according to the third and fourth 
column of Tab. respectively. Calculate [yk] and [Vy^] according to the 
third and fourth column of Tab. 1,7.71 respectively. Compute [A^j depending on 
Li, Lj, and Ck according to the second column of Tab. \A.l\ in the appendix. 

3. Compute [An-i-t] from [A|j_|_(] and Cn+t according to Eg. (14.5p and set [p] = 
[Vn+t], [V(/3] = [Vy„+t], and [A^^] = [A„+t]- 

Table 4.5 

Rules for the computation of the sets and Ck in the k-th line of the codelist lIS.lll for variables 
and binary operations (left) and compositions (right). Rules for yk are repeated here for convenience. 


op 4>A; 

Vk 

Ik 

Ck 

op 

Vk 

Ik 

Ck 

VcLT 

Xk 

JT\{k} 

N 

powNat 

vT 

li 

li 

add 

Vi + Vi 

li n ij 

Ci n Cj 

oneOver 

1/j/i 

li 

li 

mul 

ViVj 

li n Ij 

li n Ij 

sqrt 

y/¥i 

li 

li 





exp 

exp(yi) 

li 

li 





In 

ln(i/i) 

li 

li 





addC 

Vi + c 

li 

Ci 





mulByC 

cyi 

li 

Ci 


Proof. The claims (fi{x) € [(fi] and V(p{x) G [V(p] for 3 II x G B are covered 
by Thm. 13.11 It remains to prove that [Ai„in(V^i^(a:)), Aniax(V^(/ 3 (a:))] C [Ag,] for 
all X G B. Since ykix) = Xk for k = l,...,n, the functions yk(x), k G A/”, are 
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independent of Xj for every j S Ifc = A/" \ {fc} and at most linearly dependent on Xj 
for every j £ Ck = A/”. Thus, the reduced Hessian reads Vj^cykix) = y^yk{x) = 0o_o 

and [A^] = [0,0] for every k € M. Now assume eigenvalue bounds [A|], ..., [A[] for 
the reduced Hessians V^cyifa;),..., Vj-cydx) and index sets Xi,... ,Ii and Ci,... ,£[ 

1 I 

have been calculated for some I £ Nn,n+t-i, and let A: = Z + 1. Since $fe(yi, ■ • ■, 2 /fe-i) 
is one of the unary or binary functions listed in Tab. 13.11 (and therefore Tabs. lT5l and 
EtJ, it depends on either one (say yi) or two (say yi and yj) of the intermediate 
variables yi,..., yu-i- The remainder of the proof must be carried out for each type 
of operation separately. We state the proof for one of the mul cases and claim 
the remaining cases can be shown accordingly. Let g{x) = yi{x), h(x) = yj{x), and 
f{x) = yk{x), which implies f{x) = g(x)h{x), since the operation in the fc-th line 
is of type mul. In order to compute eigenvalue bounds [AJ,] for the reduced Hessian 
Vj^cf^x), we first evaluate the index sets 1/ and £/. According to Lem. 14.141 we 
obtain 


Xf =igr\ih=x^r\Xj and Cf =XgnXh =x,nXj, 

where we used Xg = Xi and Xh = Xj , which hold by construction. Assuming we have 
CgUCh C A/” and Cj — Cg C Ch, applying Lem. 14.151 fspecifically, rule 12 in Tab. l4.4p 
results in 

[A]p] = [At([V£cy],[V£=/i])] +[/i][At] +[y] [min{At,0},max{A|^,0}], 

= [At([V£=y*], [S/ciVi])] + [vj] [A1] + [yt] [min{Aj, 0}, max{A], 0}], 

where the second equation results from substituting the codelist notation [g] = [yi], 
[h] = [Vj], [Vy] = [Vy,], [V/i] = [Vy^], [Aj] = [A|], and [Aj^] = [A]]. Finally, since 
Ck = Cf and consequently V\cyk{x) — \/j-cf{x), the eigenvalues of Vj-cykix) are 

conhned to [A^j = [AJ,] for all x £ B. Since the second equation in (14.251) is equal to 
the rule in Tab. lA.ll for the case mul and Ci U £ N and Ck = C-i C Cj, this proves 
the claim for the selected case. □ 

ProDOsition l4.17l is illustrated with two examples. First, we revisit the motivating 
Exmp. [TJ Recall that we evaluated the conservative eigenvalue bounds [A,^] = [0,4] 
using the original method from [9]. 

Example 2 (improved method applied to (f{x) = xf + X 2 from Exmp. |T]). Con¬ 
sider the function cp from Exmp. [I| again. Proposition \4-i1l\ results in the following 
extended codelist. Note that we do not list the expressions for [yk] and [Vy^j in (|4.26l) 
since they are identical to the corresponding expressions in glD. Further note that 
Xk and Ck are independent of B. 


k 

Xk 


Ck 


[Al] 

1 

Af\{l} 

= {2} 

N 

= {1,2} 

[0,0] 

2 

N\{2} 

= {1} 

N 

= {1.2} 

[0,0] 

3 

Xi 

= {2} 

Xx 

= {2} 

2[A.([V£.yi[)[=2[A.([V{ijyi[)[ 

4 

X 2 

= {1} 

X 2 

= {1} 

2[A.([V£cy2[)[=2[A.([V{2}y2[)[ 

5 

Xi n 24 

= 0 

Ci n Ci 

= 0 

[min{Aj, Ai}, max{A 3 , A}}] 


[A^] = [Xl] 


The expressions for [aJ,] in lines 3 and 4 of the extended codelist in (14.261) refer to 
the first rule associated with the powNat-operation in Tab. \A.1\ since Ci = M and 
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£2 = , respectively. Since £3 = {2} C Af, £4 = {1} C Af, and £3 U £4 = Af, 


we obtain the bounds [Aj] = [min{^, A 4 }, inax{A 3 , A 4 }] according to the last rule for 
the add-operation in Tab. \A.1\. Finally, since £5 = 0, we have [A,^] = [A 5 ] = [A 5 ] 
according to Eq. (li3D . 

Evaluating the extended codelist (14.261) for the hyperrectangle B = [0, 1] x [0, 1] 
(as in Exmp. m by computing and [Vj/fc] according to (|4.1I) and [Afc] according 
to (14.2611 yields [Aj] = [A 2 ] = [0,0] and [A 3 ] = [A 4 ] = [A 5 ] = [2,2], where we used 
[A4[V{i}?/i])] = [As([l,l])] = [1,1] and [A 4 [V{ 2 } 2 / 2 ])] = [As([l,l])] = [ 1 , 1 ] (see 
Eq. (14.21) for numerical results on [yk] and [Vy^]/ Thus, using the improved method, 
we obtain the tight eigenvalue bounds [A,^] = [Ag] = [A*] = [2,2]. 

We analyze another example to demonstrate that the new method results in 
considerable improvements for all functions that involve multiplications. In fact, we 
know from [SJ Rem. 4.3] that 0 € [Ai^] for the original method if the mul-operation is 
required in the codelist of any ip with n > 2. This is a severe drawback of the original 
method, since it implies that any convex (resp. concave) function ip : R" ^ R 
involving mul-operations will never be identified to be convex (resp. concave) using 
the method from [HI . The following example shows that this restriction does not apply 
for the improved method. 

Example 3 (comparison of [5] and improved method for ip(x) = xl + X 2 exp(a; 2 )). 
Consider the function (/? : R^ —>■ R with ipix) = x\ + a;2exp(a;2) on a S C R^. 
Theorem \3.1\ (i.e., the original method from m) results in the following extended 
codelist, where the expressions for yk are only listed for illustration of the codelist dSU 
of ip. We skip the first three lines, since they are identical to those in 611). 


k 

yk 

[Vk] 

Nvk] 

[Xk] 

4 

exp(i/2) 

[exp([2/2])] 

[2/4] [Vl/2] 

[l/4]([W([Vl/2])] + [A2]) 

5 

2/21/4 

[2/2] [2/4] 

[2/4][Vl/2] + [2/2][Vj/4] 

[2/4] [A2] + [2/2] [A4] + [At ( [Vi/2], [V1/4])] 

6 

2/3 + 2/5 

[2/3] + [2/5] 

[V1/3] + [Vy5] 

[A3] + [A5] 


= 2/6 

[t] = [ye] 

[V(^] = [Vye] 

[A(p] = [Ae] 


Evaluating this codelist for B = [0,1] x [0,1] yields 
[A,^] = [Ae] = [- exp(l) + 1,3 exp(l) + 2] 


-1.7183,10.1548]. 


Proposition [£73 (i.e., the improved method) results in the following extended 
codelist. The first three lines are identical to those in (14.261) in this case. Eurther 
note that the expressions for [yk] and [Vj/fe] can be found in Eq. (|4.1I) (lines T3) and 
Eq. (14.271) (lines 4-6). 


k 

Tk Ck 

[Al] 

4 

X2 ={1} l2 ={1} 

[ 2 / 4 ] [As([V£c 1 / 2 ])] = [ 2 / 4 ] [As([V{ 2 } 2 / 2 ])] 

5 

X 2 0 X 4 = {1} X 2 0 X 4 = {1} 

[At([V£g2/2], [V£g2/4])] + [2/4] [A^] + [2/2] [aJ] 

= [At([V{2}2/2], [V{2}2/4])] + [2/4] [A^] + [2/2] [A|] 
[min{ A^ , A^ }, max {A3, Ag }] 

6 

Is n 2:5 = 0 £3 n £5 = 0 



[A^j] = [Ag] 


Evaluating (14.281) for S = [0,1] x [0,1] yields 

[A^] = [At] = [2,3exp(l)]«. [2,8.1548], 

Just as in Exmp. d the improved method results in tight spectral bounds while the 
original method from m provides loose outer approximations. In particular, 0 G 
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[—1.7183,10.1548] for the original method as predicted by {PI Rem. 4-3] but 0 ^ 
[2,8.1548] for the improved method presented here. Convexity of ip on B can therefore 
be established with the improved but not with the original method. 

More generally, the improved method results in eigenvalue bounds that are always 
as tight as, or tighter than, the original method from [9], as stated in the following 
proposition. 

Proposition 4.18 (accuracy of the improved method). Assume ip is twice 
continuously differentiable on U and can be written as a codelist (13.11) . Let B — 
[xi] X • • • X [xn] ClL be arbitrary and denote the eigenvalue bounds for 'S/‘^p(x) on B 
computed according to Thm. \4-17\ and Prop. \4-17\ by and 

respectively. Then, 



Proof. The proof immediately follows from Lems. 14.8114.111 and 14.151 □ 

In [9l Prop. 4.4] it was shown that the numerical complexity for evaluating 
the extended codelist resulting from Thm. 13.11 is of order 0{n) N{p), where N{p) 
denotes the number of operations needed to evaluate p & point in its domain. It is 
remarkable that this order of complexity can be maintained for the improved method. 
This is summarized in the following proposition. 

Proposition 4.19 (numerical complexity of the improved method). Assume p 
is twice continuously differentiable on U and can be written as a codelist iH) with 
t = N{p) operations. Let iV([v5], [Vyi], [A,^]) denote the number of operations that are 
necessary to calculate the bounds [p] C M, [Vyi] C K”, and [A(^] C K for a given 
hyperrectangle B CU using the extended codelist from Prop. |^.17[ Then, 

N{[p\,[^p\,[X,,\) = 0{n)N{p). 


Since the proof of ProD. l4.19l is very similar to that of [9l Prop. 4.4], we only sketch 
it. The extended codelist that results from Prop. |4T3 involves the index sets Ik and 
Ck, which were not required in the original method. These index sets do not depend 
on the particular hyperrectangle B as pointed out in Exmp. [2 but they are uniquely 
determined by the function p itself. Consequently, all index sets need to be determined 
only once. This step can be carried out at the time of construction of the extended 
codelist. In particular, it need not be repeated at the time of evaluating the codelist 
for a particular B. Once Ik and Ck have been determined, each line of the extended 
codelist is specified by the rules in Tab. lA.ll fand Tab. 13.11) . It is easy to show that the 
evaluation of every expression in the second column of Tab. lA.U rea uires at most 0{n) 
basic operations (like additions, multiplications, or comparisons of two real numbers; 
see [9l Sect. 4.1] for further details). Thus, under the assumption that [(^] and [Vf/j] 
are known, we need 0{n)N{p) basic operations for the computation of [A,^]. Since 
the calculation of [p] and [Vtp] require 0(1) and 0{n) basic operations according 
to standard results in AD and lA (see, e.g., 0 0), we obtain lV([</3], [Vtp], [A,^]) = 
0(1) N{p) + 0(n) N{p) + 0(n) N{p) = 0{n) N{p). 

5. Numerical experiments for a large number of examples. In this sec¬ 
tion, we analyze 1522 numerical examples taken from the COCONUT collection of 
optimization problems [T3]. We consider all COCONUT problems with 1 < n < 10 
variables and extract those cost and constraint functions that can be decomposed into 
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the operations listed in Tabs. IXT] and lA.fl For each function ip : R." —>■ K, we consider 
100 (randomly generated) hyperrectangles 6 C in the domain T) oi p specified in 
the respective COCONUT problem. For ease of comparison, the set of examples as 
well as the associated hyperrectangles are identical to the examples considered in m- 

For each of the resulting 1522 • 100 sample problems, we solve problem (EB 
using the improved algorithm (A^ for short) in Prop. 14.171 We compare the resulting 
eigenvalue bounds with those obtained from two established methods using interval 
Hessians (see problem (11.31) 1 and either Gershgorin’s circle criterion (G for short) 
or Hertz and Rohn’s method (H for short) for the computation of spectral bounds of 
interval matrices (see problem (11.41) 1. We choose G and H as reference procedures due 
to the favorable computational complexity of G and since H provides tight eigenvalue 
bounds for problem (11.41) (cf. Sect.[T|). We refer to the original papers UlSllllIl] or 
the summaries in [HI HJ] for a detailed description of methods G and H. 

Table 5.1 

Classes used to aggregate results in Tab. 15.21 Symbols [A^tl = [-^g] = [Ag’-^g], md 

[A/f] = [Aifi denote eigenvalue bounds calculated by the improved algorithm A} (see Prop, \4-17\ l, 
Gershgorin’s circle criterion and Hertz and Rohn’s method, respectively. 


class 

verbal definition 

formal dehnition 



lower bound A 

upper bound A 

1 

A^ worse than G (and H) 

AaI < Ag ^ Ah 

Ah < Ag < A^t 

2 

A^ equal to G but worse than H 

Ag ~ AaI < Ah 

Ah < A^t ~ Ag 

3 

A^ better than G but worse than H 

Ag < AaI < Ah 

Ah < A^t < Ag 

4 

A^ equal to H (and equal to or better than G) 

Ag ^ Ah ~ AaI 

A^t ~ Ah < Ag 

5 

A^ better than H (and G) 

Ag ^ Ah < AaI 

Aai < Ah < Ag 


For each sample problem, we analyze whether A^ performs better than, equally 
good as, or worse than G and H. We independently compare the lower and upper 
eigenvalue bounds of the particular methods and categorize the results according to 
the five classes in Tab. EH Note that G never performs better than H (since H 
provides tight bounds for (11.41) 1. Consequently, the relations Aq < Ajj and Ah < Ag 
always hold. Hence, the list of classes in Tab. I5.l| is complete in the sense that every 
example can be uniquely classified into one of the five classes. It remains to comment 
on the precise meaning oi a > b and a ~ 6 as used for the classification in Tab. 15.11 
To this end, we introduce the function 


dev(a, b) 


a — b 

1 + 0.5 \a + b\ 


which evaluates a weighted difference oi a, b G R. Based on dev(a, b), we define 


a > b 


dev(a, b) > e and 


|dev(a, 6)1 < 


(5.1) 


where e G R+ represents an error bound. Note that |dev(a, 6 )| is approximately equal 
to the relative error for two large but almost equal numbers a, 6 G R and almost equal 
to the absolute error for two small but almost equal numbers a,b gM.. This behavior 
is useful since the absolute values of the computed eigenvalue bounds range across 
multiple magnitudes. 

We summarize numerical results for the analyzed examples in Tab. 15.21 (with 
e = 10“®). We list the percentage of samples that fall into the classes 1 to 5 from 
Tab. 15.11 separated by dimension n of the underlying example. In order to compare 
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Table 5.2 

Numerical results for 1522 • 100 sample problems. The classification of the results is carried out 
according to Tab. \ 5. f\ and Eq. (15.11 1 with e = 10“®. For each class, the left respectively right column 
refers to eigenvalue bounds obtained by the original algorithm A from m and the improved variant 
taking sparsity into account. 


examples 
n # 

1 

A 

Al 

2 

A 

At 

3 

A 

At 

A 

4 

At 

5 

A 

At 

2 

62 

56.37 

29.06 

1.07 

1.08 

12.24 

16.06 

18.89 

41.27 

11.43 

12.54 

3 

1078 

2.98 

0.95 

77.53 

0.31 

0.85 

1.01 

17.70 

96.72 

0.95 

1.01 

4 

67 

60.75 

35.28 

4.84 

4.63 

8.45 

16.76 

15.04 

32.14 

10.92 

11.19 

5 

88 

56.80 

35.34 

3.32 

0.21 

10.16 

14.90 

14.91 

34.47 

14.81 

15.09 

6 

95 

35.05 

25.65 

5.31 

4.13 

31.91 

34.25 

13.87 

21.38 

13.86 

14.58 

7 

27 

65.80 

27.07 

11.93 

8.50 

0.02 

19.44 

22.26 

44.98 

0.00 

0.00 

8 

15 

94.23 

63.83 

3.20 

4.23 

1.27 

29.80 

1.30 

1.53 

0.00 

0.60 

9 

24 

57.27 

25.02 

8.29 

0.02 

18.71 

34.94 

4.25 

22.58 

11.48 

17.44 

10 

66 

51.17 

12.24 

1.00 

0.22 

22.12 

34.11 

15.86 

41.54 

9.86 

11.89 

all 

1522 

17.77 

9.09 

56.11 

0.94 

5.32 

7.78 

16.86 

77.90 

3.95 

4.30 


the improved algorithm in Prop. ITTTi to the original method from [9] (see Thm. IXTll . 
we also list the classification results using the original algorithm (A for short). The 
numerical results confirm that the consideration of sparsity significantly improves the 
tightness of the computed eigenvalue bounds. To see this, note that for each dimension 
n, the percentages in class 1 (where the established approaches outperform the direct 
computation of eigenvalue bounds) decrease while the percentages in classes 4 and 
5 (where the direct computation of eigenvalue bounds performs as good as or better 
than Hertz and Rohn’s method) increases using the improved algorithm instead of 
the original A. In particular, it is remarkable that the improved algorithm A^ results 
in worse eigenvalue bounds than G in only 9.09% of all cases in contrast to 17.77% 
for the original method A. Moreover, A^ provides equally good or better eigenvalue 
bounds than H in 82.20% = 77.90% + 4.30% of all cases while the corresponding 
percentage only reads 20.83 % = 16.86 % + 3.95 % for A. 

Another observation is that the ratios in the particular classes seem to be indepen¬ 
dent of the dimension n (i.e., there is no trend). This is important since the numerical 
complexities of the established approaches G and H vary between 0{n) N{ip) -|- 0{'n?) 
and 0{'n?) N{ip) -|- 0(2" n^) operations (see Sect. [Hand the benchmark in [E]), while 
the direct eigenvalue bound computation requires 0(n) N(ip). Thus, methods A and 
Al become numerically very attractive for high dimensions n. 

Table 5.3 

The last line of Tab. \5.^ for different choices of the error bound e in ISTl l. 


e 

1 

A 

At 

2 

A 

At 

3 

A At 

4 

A At 

5 

A At 

10 "^ 

17.16 

8.78 

54.11 

0.91 

4.29 

6.57 

20.78 

79.74 

3.66 

4.01 

10 -® 

17.77 

9.09 

56.11 

0.94 

5.32 

7.78 

16.86 

77.90 

3.95 

4.30 

10-7 

18.04 

9.28 

57.62 

0.66 

5.61 

8.24 

14.67 

77.41 

4.06 

4.41 


According to (EB, the classification in Tab. 15.21 depends on the choice of the 
error bound e. We repeated all calculations for various choices of e and present the 
results reported in the last line of Tab. 15.21 for e = 10“® and e = 10“^ in Tab. 15.31 
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As expected, the ratios in classes 1 and 5 increase for decreasing e, since we detect 
Aaj < Ag < ^Ati Ah < AAt) and A^t < Ah) for a larger number of 

examples (cf. (15.11) 1. However, beside this observation, the results are robust w.r.t. 
the value of e. 

6. Conclusion. We significantly improved a method recently introduced in [9] 
for the efficient computation of spectral bounds for Hessian matrices of twice contin¬ 
uously differentiable functions on hyperrectangles. The improvements build on the 
identification and utilization of sparsity that naturally arises in the first lines of every 
codelist for a function tf : K." —>■ R. 

The improved method was applied to a set of 1522 examples previously analyzed 
in m- The numerical results show that the consideration of sparsity results in signifi¬ 
cantly tighter eigenvalue bounds. In fact, the improved method provided equally good 
or better eigenvalue bounds than Hertz and Rohn’s method in 82.20 % of the examples 
while the corresponding percentage only reads 20.83 % for the original procedure. 

In addition to illustrating the practical usefulness of the proposed improvements, 
we provided an important theoretic result. In fact, it is well-known that the original 
method from [9] results in spectral bounds with 0 € [A,^] for any function that involves 
the multiplications of two or more variables (see O Rem. 4.3]). Consequently, convex 
functions that involve such a multiplication cannot be detected to be convex with the 
original method. We showed that this restrictions does not apply for the improved 
method. 

Acknowledgements. Funding by Deutsche Forschungsgemeinschaft grant MO- 
1086/9 is gratefully acknowledged. 
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Appendix A. Supplementary results. 

Lemma A.l. Let [a], [ 6 ], [c] C M. be real intervals, let [A*([a], [ 6 ], [c])] be defined 
as in (14.211) . and consider the matrix set 

T~L = {H G I Hii s [a], H 22 S [ 6 ], H 12 G [c], H = 


Then, 


A*([a],[6],[c]) = mjn An 




and 


max Amax 
H£H 


{H)=A.{[a],[b],[c]). 


(A.l) 


Proof. The eigenvalue bounds of a symmetric matrix H = {%l) read 


— - 

(^a + b- y/{a- 

bY 

+ 4c2j , 

'^max(-^) — 2 

+ 6 + ■\/(a — 

bY + 

We therefore 

have to show that 





a + b — 

Via-kY + d = 


min 

a + b— — 

&)2 + 4c 2 and 

(A.2) 


aG[ 

a\,bG[b\,cG[c\ 

a b 

\J {a — bY + d = 


max 

a + b+ Y{a — 

bY + 4c2, 

(A.3) 


oG[ 

a\,bG[b\,c£[c\ 


where the l.h.s. results from (I4.21|) and where d = A max{c^, c^}. We show that (IA.2|) 
holds and claim (|A.3I) can be proven analogously. First note that the r.h.s. in (jA.2|) 
can be simplified to 

min a-\-b— (a — bfi + 4c^ = min a + b — (a — b)"^ + d (A.4) 

aG[a],hG[b],cG[c] aG[a].bG[6] 

since c only occurs in the radicand. Consider the function / : ^ K, f{a,b) = a + 

b— \/{a — 6)2 + d, which occurs on the r.h.s. of (IA.4I) and note that /(a, b) is concave 
(since g(a,b) = ^J{a — bY -\- d is convex). Since the hyperrectangle B = [a\ x [ 6 ] is 
convex, the minimum on the r.h.s. of (IA.4I) is attained at one of the vertices of B. 
Among the candidate tuples (a, 6), {a,b), {a,b), and (a, &), it is easy to show that 
(a, 6 ) results in the smallest function value, i.e., f{a,b) < min{/(a, &), f{a,b), f{a,b)}. 
Thus, (|A.2I) holds. □ 
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Table A.l 

Rules for the calculation of the eigenvalue bounds [A^] in the k-th line of the codelist O.lll 
according to Prop. |4 -?7j The intervals [At], [At^o] shorthand notations for [At] = 

[At([V£cyi], [V£cyj])], [Ai,o] = [min{A|,0},max{A|,0}], and[Xjfi] = [min{A], 0}, max{A], 0}]. The 
index sets Cu and Cn are shorthand notations for Cu = CiU Cj and Cn = CiH Cj. Condition C* 
reads {li U Ij = M) A {\Ii\ = n — 1) A {\Ij | = n — 1). 


op 

[All 

condition 

add 

[0,0] 

Ci=M hCj^M 


[All 

Ci (Z N h Cj = M 


[A]] 

Ci=NACj zN 


[A.([Al],[At])] 

Ci z N A Cj z N A Cu = N 


[a1] + [aJ] 

CuZNA Ci= Cj 


[a1] + [Aj,o] 

Cu ZN ACiZ Cj 


[Ai,o] + [a]] 

CuZN A Cj Z Ci 


[Ai,o] + [Aj,o] 

Cu zN ACi <^Cj ACj ^Ci 

mul 

[At] 

Ci=M ACj=M 


[At] + [y,] [All 

Ci z N A Cj = N A Ck = Ci 


[At] + [Vj] [Ai,o] 

Ci z N A Cj ^ N A Ck z Ci A -.C* 


[A*([yj] [A|], [0,0], [Vici/i][Vict/j])] 

Ci zNACj = A/'A C* 


[At] + [yi] [A+] 

Ci =NACj zNACk = Cj 


[At] + [yi] [Aj,o] 

Ci = M A Cj Z M A Ck Z Cj A -.C* 


[A*([0,0],]j/,][A]],[Vicy,][Vi.y,-])] 

Ci=M ACj c a/'ac* 


M + [A.([A.([%] [All,]?/,] [A]])]. [0,0])] 

Ci Z N A Cj Z N A Cu = N A Ck Z Cn 


[At] + [K{[y,] [Aj], [yi] [A]])] 

Ci Z N A Cj Z N A Cu = N A Ck = Cn A ^C^, 


[A*(bj] [A1], [yi] [A]], [Vicyi][Vxcy^-])] 

CiZM ACj zM AC., 


[At] + M [All + [y.] [A]] 

Cu zMACk =Ci=Cj 


[At] + [yj] [At] + [yi] [Aj.o] 

CuZNACk = CiZ Cj 


[At] + [yj] [Ai,o] + [yi] [A]] 

Cu z N A Ck = Cj z Ci 


[At] + [A.([y,] [A|] + M[A]],[0,0])] 

CuZNACk zCi = Cj 


]At] + [^r{[yj] [A|] + [yi] [Aj,o], [0,0])] 

Cu zN ACk Z CiZ Cj 


[At] + [Ar([yj] ]Ai,o] + [yi] [A]], [0,0])] 

Cu Z N A Ck Z Cj Z Ci 


[At] + [yj] [Ai,o] + [yi] [A^^o] 

Cu zM ACi '^CjACj <^Ci 

powNat 

m{m-l)[y,r-^[Ks{[Vciyi])] 

Ci=Af 


rn[yi]^~^{{m-l)[Ks{[Vciyi])] + [yi] [Aj]) 

CiZN A Ck=Ci 


™-[2/i]'""^((”i’-l)[As([V£ji/i])] + [yi] [Ai,o]) 

Ci zAf ACk ZCi 

oneOver 

^[ykY [As([V£=i/i])] 

Ci=JV 


M^(2[y,][A.([V£ct/,])]-[Al]) 

Ci zAf ACk ^ Ci 


[VkY (2 [t/fc] [As([V£c j/i])] - [Ai,o]) 

Ci Z Af A Ck Z Ci 

sqrt 

l/(-4[t/fe]^)[A4[V£=t/t])] 

Ci=Ar 


1/(2 ]j/fe])(l/(-2[t/,])[A.([V£=t/.])] + [Al]) 

CiZAT ACk =Ci 


1/(2 [ 2 /fc])(l /(-2 [2/,])[A,([V£=t/,])] + [At,o]) 

CiZAf ACk zCi 

exp 

[Vk] [As([V£=2/i])] 

Ci=Af 


[Vk] ([A.([V£=yt])] + [Al]) 

CiZAf ACk =C. 


[Vk] ([Aa([V£gyi])] + [Ai,o]) 

CiZAf ACk zCi 

In 

“l/[j/i] [As([V£= 1/i])] 

Ci=Af 


l/N([Al]-l/N[A.([V£=t/.])]) 

CiZN ACk = Ci 


l/[j/i] ([Ai,o] - l/[yi] [As([V£=yi])]) 

Ci Z Af A Ck Z Ci 

addC 

[0,0] 

Ci=Af 


[Al] 

Ci zAf 

mulByC 

[0,0] 

Ci=Af 


c[A|] 

CiZAf 
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